cap program drop RDplot
program define RDplot

syntax varlist(min=2 numeric) [if], [wid(real 0) h(real 0) z0(real 0) z1(real 0) Degree(integer 0) Kernel(string) BW(real 0) PW(real 0) Ylab(string) Xlab(string) Ytitle(string) Yformat(string) save(string)]

//Default parameters
	local y_format "%9.2f"
	if "`yformat'" != "" {
	local y_format `yformat'
	}


//Setup
	tempname yvar run_var running bins seq means x Ly_hat Lse Ry_hat Rse y_hat se CI_low CI_up
	qui gen `yvar' = `1'
	qui gen `run_var' = `2'
	
//Evenly Spaced Bins
	qui gen `running' = .
	qui gen `means' = .
	
	local j = 1
	qui foreach i of numlist `z0'(`h')`z1' {
	
		//Width
		if `i' < 0 {
			local w0 = `i' - `wid'
			local w1 = `i' + `wid'		
			if `w1' >= 0 {
				local w1 = 0
			}
		}
		if `i' >= 0 {
			local w0 = `i' - `wid'
			local w1 = `i' + `wid'
			if `w0' < 0 {			
				local w0 = 0
			}
		}
		
		//Local Means
		qui sum `yvar' if `run_var' >= `w0' & `run_var' < `w1'
		local mu = round(r(mean), 0.0001)
		local N = r(N)
			
		replace `running' = `i' in `j'
		replace `means' = `mu' in `j'

		*dis in yellow "`i' (Obs = `N') [`w0' ,`w1'): `mu'"		
		
		local j = `j' + 1
	}
	//qui replace `running' = `running' + `h'/2 if `running' == 0
	
//Local Polynomial Regression
	//Grid
	qui gen `x' = .
	local i = 1
	qui foreach j of numlist `z0'(`h')`z1' {
	replace `x' = `j' in `i'	
	local i = `i' + 1 
	}
	
	//Regression Options
	local lpoly_ops degree(`degree') kernel(`kernel') bw(`bw') pw(`pw') n(10) nograph
	
	//Regressions at both sides of the cutoff
	lpoly `yvar' `run_var' if `run_var' < 0, generate(`Ly_hat') se(`Lse') at(`x') `lpoly_ops'
	lpoly `yvar' `run_var' if `run_var' >= 0, generate(`Ry_hat') se(`Rse') at(`x') `lpoly_ops'
	qui gen double `y_hat' = `Ly_hat' if `x' < 0
	qui replace `y_hat' = `Ry_hat' if `x' >= 0
	
	//Confidence Intervals
	qui gen double `se' = `Lse' if `x' < 0
	qui replace `se' = `Rse' if `x' >= 0
	qui gen `CI_low' = `y_hat' - abs(invnorm(0.025))*`se'	
	qui gen `CI_up' = `y_hat' + abs(invnorm(0.025))*`se'	
	
//Plot 
	# delimit ;
	tw  (rarea `CI_low' `CI_up' `x' if `x' < 0, fcolor(gs13%40) alw(none) clcolor(black) clwidth(medthin))
		(line `y_hat' `x' if `x' < 0, lcolor("0 0 200"))
		(rarea `CI_low' `CI_up' `x' if `x' >= 0, fcolor(gs13%40) alw(none) clcolor(black) clwidth(medthin))
		(line `y_hat' `x' if `x' >= 0, lcolor("0 0 200"))
		(scatter `means' `running', mfcolor(gs10) mlcolor(black) mlw(thin) msize(medsmall))
		,
		graphregion(color(white)) plotregion(lcolor(black))
		ylab(`ylab', format(`y_format') labsize(medsmall) angle(horizontal) grid gmin gmax glwidth(thin) glcolor(gs10) glpattern(dot) nogextend) 
		xlab(`xlab', format(%9.1fc) labsize(small) grid gmin gmax glwidth(thin) glcolor(gs10) glpattern(dot) nogextend)
		xti(" " "Distance to National Distinction Award cutoff ({&sigma})") 
		yti("`ytitle'" " ") legend(off)
	;# delimit cr
	
	if "`save'" != "" {
	gr export "`save'", replace wid(2000) hei(1500)
	}
	
end

